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Abstract 

The multi-level Monte Carlo method proposed by M. Giles (2008) approximates the 
expectation of some functionals applied to a stochastic process with optimal order of 
convergence for the mean-square error. In this paper, a modified multi-level Monte 
Carlo estimator is proposed with significantly reduced computational costs. As the 
main result, it is proved that the modified estimator reduces the computational costs 
asymptotically by a factor (p/a) 2 if weak approximation methods of orders a and 
p are applied in case of computational costs growing with same order as variances 
decay. 
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1 Introduction 



The multi-level Monte Carlo method proposed in [7] approximates the ex- 
pectation of some functional applied to some stochastic processes like e. g. 
solutions of stochastic differential equations (SDEs) at a lower computational 
complexity than classical Monte Carlo simulation, see also [5,8,9]. Multi-level 
Monte Carlo approximation is applied in many fields like mathematical finance 
[1,6], for SDEs driven by a Levy process [3] or by fractional Brownian motion 
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[11]. The main idea of this article is to reduce the computational costs fur- 
ther by applying the multi-level Monte Carlo method as a variance reduction 
technique for some higher order weak approximation method. As a result, the 
computational effort can be significantly reduced while the optimal order of 
convergence for the root mean-square error is preserved. 

The outline of this paper is as follows. We give a brief introduction to the 
main ideas and results of the multi- level Monte Carlo method in Section 2. 
Based on these results, in Section 3 we present as the main result a modified 
multi-level Monte Carlo algorithm that allows to reduce the computational 
costs significantly. Depending on the relationship between the orders of vari- 
ance reduction and of the growth of the costs, there exists a reduction of the 
computational costs by a factor depending on the weak order of the underly- 
ing numerical method. As an example, the modified multi-level Monte Carlo 
algorithm is applied to the problem of weak approximation for stochastic dif- 
ferential equations driven by Brownian motion in Section 4. 



2 Multi-level Monte Carlo simulation 

Let (Q, J 7 , P) be a probability space with some filtration (J : t)t>o and let X = 
(Xt)tei denote a stochastic process on the interval / = [to,T] adapted to the 
filtration. In the following, we are interested in the approximation of Ep(/(X)) 
for some functional / € F where F denotes a suitable class of functionals that 
are of interest. Further, let an equidistant discretization I h = {t ,ti, . . . ,t^} 
with < t < tx < . . . < t N = T of the time interval I with step size h 
be given. Then, we consider a probability space (Cl, J 7 , P) with some filtration 
{J-t)tei h an d we denote by Y = (Y t ) te i h a discrete time approximation of X 
on the grid Ih, adapted to (jF t )tei h - Here, the probability spaces (f^J 7 , P) 
and (fl, J 7 , P) may be but do not have to be equal and we assume that Y 
approximates X in the weak sense with some order p > 0, i.e. 

||Ep(/(F))-E P (/(X))||=0(^) (1) 

for all / e F. 

In order to approximate the expectation of f{X) we apply the multi-level 
Monte Carlo estimator introduced in [7]. For some fixed MeN with M > 2 
and some L e N we define the step sizes hi = and let Y l = (Y t ) te j h denote 
the discrete time approximation process on the grid 1^ based on step size hi 
for I — 0, 1, . . . , L. Then, the multi- level Monte Carlo estimator is defined by 

Yml = EY 1 (2) 

1=0 
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for some L G N using the estimators Y° = ^- /(^ oW ) and 
for I = 1, . . . , L. Then, we get 

e p (y ml ) = Ep(/(r )) + EEp(/(y') - Z^- 1 )) • (4) 

Z=l 

Here, we have to point out, that both approximations Y 1 ^ and Y l ~ 1 ^ are sim- 
ulated simultaneously based on the same realisation of the underlying driving 
random process whereas (Y l ^\ Y l ~ 1 ^) and (Y l ^\ Y 1-1 ^) are independent re- 
alisations for % ^ j. 

Now, there are two sources of errors for the approximation. On the one hand, 
we have a systematical error due to the discrete time approximation Y h based 
on step size h which is given by the bias of the method. On the other hand, 
there is a statistical error from the estimator for the expectation of f(Y) by the 
Monte Carlo simulation. Therefore, we consider the root mean-square error 

e(Y ML ) = (E P (\Y ML - Ep(/(X))| 2 )) V2 (5) 

of the multi-level Monte Carlo method in the following. In order to rate the 
performance of an approximation method, we will analyse the root mean- 
square error of the method compared to the computational costs. Therefore, 
we denote by C(Y) the computational costs of the approximation method Y. 
In order to determine C(Y), one may use a cost model where e.g. each oper- 
ation or evaluation of some function is charged with the price of one unit, i.e. 
one counts the number of needed mathematical operations or function evalu- 
ations. Further, each random number that has to be generated to compute Y 
may also be charged with the price of one unit. 

Then, it is well known that the optimal order of convergence in the case 
of the classical Monte Carlo estimator Y MC = j^Yh =1 f(Y^) is given by 
g{Ymc) — O {{1/C{Ymc)) 2p+1 ) where p is the weak order of convergence of 
the approximations Y, see Duffie and Glynn [4]. Thus, higher order weak ap- 
proximation methods result in a higher order of convergence with respect to 
the root mean-square error. Clearly, the best root mean-square order of con- 
vergence that can be achieved is at most 1/2. However, the order bound 1/2 
can not be reached by any weak order p approximation method in the case of 
the classical Monte Carlo simulation. Therefore, in order to attain the optimal 
order of convergence for the root mean-square error we apply the multi-level 
Monte Carlo estimator (2). The following theorem due to Giles [7] is presented 
in a slightly generalized version suitable for our considerations. 
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Theorem 2.1 For some L e N 7 /ei denote the approximation process on 
the grid 1^ with respect to step size hi — for each I — 0, 1, . . . , L, respec- 
tively. Suppose that there exist some constants a > and Ci jQ , c 2j o, c 2 , c 2) l > 
and (3, (3l > such that for the bias 

1) |E P (/(X))-Ep(/(^))|<c 1 , a / i 2 
and for the variances 

2) Varp(/(y°(T)))<c 2i0 ^ ; 

3) Var p (/(r i (T)) - f(Y l ~ 1 (T))) <c 2 h? for I — 1, . . . , L — 1, 

4) Varp(/(y i (T)) - f(Y L ~\T))) < c 2 , L h{\ 

Further, assume that there exist constants c 3t0 ,c 3 ,c 3j L > and 7,7l > 1 such 
that for the computational costs 

5) C(Y°(T))<c 3fi Th^, 

6) C(Y l (T), Y l -\T)) < c 3 T V for I — 1, . . . ,L — I, 

7) C(Y L (T),Y L -\T)) < c 3 , L Th L ^. 

Then, for some arbitrarily prescribed error bound e > there exist values L 
and Ni for I = 0, 1, . . . , L, such that the root mean-square error of the multi- 
level Monte Carlo estimator Y ML has the bound 

e(Y ML ) < e (6) 

with computational costs bounded by 

c 4 e~ 2 if (5 > 7, p L > 7 L , a > \ max{7,7 i }, 

C(Y ML ) < I c 4 e~ 2 (log(e)) 2 if (3 = 7, (3 L > >y L , a > \ max{7,7 L }, 

c 4 £- 2 - maJC{7 "^" fa} if /3 < 7, a > Wml-^-fa-fel 

(7) 

for some positive constant c 4 . 

The calculations for the proof follow the lines of the original proof due to Giles 
[7]. Considering the mean square-error 

e(Y ML ) = (I Ep(/(X)) - E P (/(F L ))| 2 + Var P (F ML )) 1/2 < e (8) 

we make use of the weight q G ]0, 1[ and claim that 

I E P (/pO) - E p (/(r L ))| 2 < qe 2 and Var P (Y ML ) < (1 - q) e 2 . (9) 

Then, we can calculate L from the bias and we have to solve the minimization 
problem 

min C(Y ML ) (10) 
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under the constraint that Vaxp(Y ML ) < (1 — q) e 2 . As a result of this, we obtain 
the following values for L and Nf. 



Let 



and N 



logiq-^c^e- 1 ^) 



a log(M) 



/3+7 / x I 
1-9 U \C3,0/ 



for I = 1,...,L - 1 and N L 
with: 



1 -2,^/C 2 ^ 



//; [ — | /■ 



1-9 

1-9 ^ \C3,i/ 



(12) 

for some g e ]0, 1[ 



In case of (3 > 7 and /3 L > j L or in case of (3 < 7 and 7^ — (3 L < 7 — j3: 



k = (c 2 , c 3 , ) 2 T 2 + (c 2 c 3 ) 2 



• In case of j3 — 7 and /?l > 7^: 



1 fa , ^aiM-^V-V 



0L-JL 



1 - Mt 



+ (c 2 ,lc 3 ,l)^ l 2 • (13) 



1 11 ^L-TL 

K = (C2, C 3 ,0) 5 + (£ - 1) (C 2 C 3 ) 5 + (c 2jL C 3iL )2 /i £ 2 



(14) 



3 The improved multi-level Monte Carlo estimator 



The order of convergence of the multi- level Monte Carlo estimator Yml given 
in (2) is optimal in the given framework. However, the computational costs 
can be reduced if a modified estimator is applied. As yet, the estimator Yml 
is based on some weak order a approximations Y l for I = 0, 1, . . . , L on each 
level. Now, let us apply some cheap low order weak approximation Y l on levels 
/ = 0, 1, . . . , L — 1 combined with some probably expansive high order weak 
approximation Y L on the finest level L. The idea is, that the approximations 
Y l contribute a variance reduction while the approximation Y L results in a 
small bias of the multi-level Monte Carlo estimator, thus reducing the number 
of levels needed to attain a prescribed accuracy. 

Let Y be an order a weak approximation method and let Y be an order 
p weak approximation method applied on the finest level. Further, let L = L p 
with 

log^d,^- 1 ^) 



L p — 



p log(M) 



(15) 
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denote the number of levels in order to indicate the dependence on the weak 
order p. Then, we define the modified multi-level Monte Carlo estimator by 



Y ML(a , p) = Y,Y l (16) 



=o 



with the estimators Y l for I — 0, 1, . . . , L p — 1 based on the order a weak 
approximations Y l as defined in Section 2, however now applying the modified 
estimator 



1 N Lp 

y Lp = -Tj— (f(Y L » f } - f{Y L ^) (17) 

Lp 1=1 

which combines the weak order a approximations Y Lp ~ l with the weak order p 
approximations Y Lp . Clearly, all conditions of Theorem 2.1 have to be fulfilled 
for Y L replaced by Y L . Then, in the case of p > a, the improved multi- level 
Monte Carlo estimator YML(a, P ) features significantly reduced computational 
costs compared to the originally proposed estimator Y ML = Y ML ( a>a y 

Proposition 3.1 Let conditions l)-7) of Theorem 2.1 be fulfilled and suppose 
that there exist constants 63,0,63,63^,^ > and 630,63^,63^ > such that 
for the computational costs 

5') C(Y°(T)) = c 3 ,o T V + ELi 4 T ho ^ , 

6') C(Y l (T), Y l ~\T)) = c 3 T hp + £* =1 4° T ^ 7+5i forl — l,...,L p — l, 
T ) C{Y^{T),Y Lp -\T)) = c 3 , Lp Thl~^ +Et 1 ^ Lp Th- L ^ +5 > 

with some 7,7 l p > 1 snc/i that^—bi > 1 and r )L p —5i > 1. Then, the multi-level 
Monte Carlo estimator Yml(oi,p) based on a weak order a > approximation 
scheme on levels 0, 1, . . . , L p — 1 and some weaA; order p > a approximation 
scheme on level L p has reduced computational costs: 

i) In case of j3 > 7 and (3 — 7 < /?l p — 7l p; inere exists some e > snc/i 
£/ia£ /or a// £ G ]0, £0] holds 

C(Y M L(a,a)){£) > 1 ^ 
C(Y ML{a:P) )(e) 

provided that a > | ; p > | max{7, 7 ip } and p > \ max{/3 + 7, /3 — 7 + 
27 ip }. In case of j3 > 7 and /3 — 7 = /?l p — 7^ p then (18) holds if in 

■^\2„_ _ „_ . _J A2c 2 ^ /1 ^^^^ c 2,%> 



addition c 2 c 3 > (1 - M 2 ) 2 c 2 , Lp C3,L p and c\f z > (1 - M 2 y c - L ^——. 

Further, for < f3 - 7 < /3 Lp - 7 Lp # no/ds C(y ML(aiP) )(£) = 0(e~ 2 ) if 
a > andp > | max{7, 7 Lp }. 
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ii) In case of f3 — 7 and (5l v > 7l p and if p > \ max{7, 7l p }, a > ], w ae£ 

lim c(W, a) )( g )^/py (19) 



^ C{Y ML ^ p) ){e) ~ \a 
and C(Y M L( a ,p))(£) = O (e~ 2 (\og(e)) 2 ) if a > and p > \ max{7, 7 Lp }. 
m) In case of [5 < 7 and 7 — /3 = 7l p — we obtain 

Um C(yML(p, P ))(g) > M 2( 7 -/3) / g 3C 2 + C 3 (C 2 C 3 ,L P ) 1/2 _ 1 

C(Y M L(a,p))(e) ' \C 3 ,L P C2,L P C 3 , Lp (c 2yLp C 3 y/ 2 V 

V2 . ^ 



(20) 

i/p > |(max{7, 7l p } — 7 + /3). // t/ie parameter q e]0, 1[ is chosen as 

« = ^ih < 21 > 

£/ien £ne computational costs C(Y ML / ap \) are asymptotically minimal. In 
general, if (3 < 7 or if f3 Lp < ^ Lp then it holds that C(Y ML ( a ^)(s) = 

max {7-/3, 7£ } 

C>(^ 2 ^ ~) /orp> |(max{7,7 Lp } -min{7-/3,7 Lp -/3 Lp }). 

Proof. Assume that e < 1. Let 5o = 0, c^g = £3,0, c 3 ^ = c 3 and cf^ Lp = c^l v - 
Then, the computational costs for Yml(q,p) are 

c{Y ML{a , p) ) = e ago t iVo + E a 3 } t V 7+51 ^ 

i=o i=o i=i ( 22 ) 



+ E4XT^^iv ip 



with L = L p 



i=0 
plog(M) 



and iVj for / = 0, 1, . . . , L p given in (12). 

Without loss of generality, suppose that 5i 7^ 5j for i ^ j and that 5k = 

4k) 4k) 4k) 
with c y 3 J = c 3 = c 3 £ p = 

use of the two estimates 



with Cg fc o = = c£l p = in the case of /3 > 7. In the following, we make 



log(e x ) logjg ^c lia T Q ) 

a " a log(M) + a log(M) ' 1 ' 

T _ 1 < MO 4. lo g(g^ c i,P TP ) , 94 x 

p -plog(M) + plog(M) • 1 } 
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Let (3 7^ 7. Then, we obtain the lower bound 



Tks 



-2 k 



c7(rML ( «,a))(e) > E ^0 2 eg, ^ + EV ' ( - ) 



> 



T 



_-2 



1-9 



1/2^ 



i=0 

fe-1 



+ ge? ( C2C2 '° C3 '° ) 1/2 t^ 

;=n V C3 / 



+ 5: 



(fe) /^ C 2 C 2,0 C 3,0 ^ 



1/2 



T^ 2 



M 2 - 1 
log(e _1 ) log(g _ 5c ljQ T c 



alog(M) cdog(M) 



2-3 



+ Q c 2 



+ 



MV-1 lalog(M) alog(M) 



(25) 



where c£ } La = c 3 i} , c 2)Lc( = c 2 , c 3jLa = c 3 , /3 Lq = /3 and 7 Lq = 7 for f M L(a,a)- 
Next, we calculate for the case of (5 7^ 7 the upper bound 

1/2 



T 2 fc 

r U — n \ 



i=0 



C 3,L P 1 



/ \ 1/2 Lp-1 

c 3,0, 
v 1/2^ 



c 3,0 



+ E V 



+h L 2 



< 



+ T E 

i=0 

T 



1-9 



7+<5* _l AW ^~^P +<5t 

1/2 



, C 3,L P / 

iW/.^isB \h /.-- .") 7-" 

-3,0' i "I" c 3 Ai « ~r c 3,L p a L p 

1=1 

i=o V ^ C3 '° ' 

+ ^ 2 ,oc 2iLp c 3 , Lp y /2 T ^ + ^^^E_^p 

+ E # (^^) T V Aj + g c« C2 A lAo 
~L V c 3 / 



A T- 



+<5; 



i=0 

+ / g (,) c2Ao + 6 ( fc) ^c^c^ 



i=0 

1/2 P Lp ~Y Lp 
Ijp 
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/ C2,oc 3fi c 2 \ 1/2 rr i^>\ / log(e x ) , log(g 2 C i iP T p ) 



V c 3 / / ^plog(M) plog(M) 

j=0 ° 3 



+ E 41 ( C2 '°^f' Lp ) 1/2 t^k^* 



=0 



i=0 \ \ C 3,L P I 

J 2^ | C 3,0 J + c 3 J _ j^-fc h C 3,L p n L p 



i=0 



(26) 



# — -y /3 — 7 | c 

(M" 1 T) i V 1 + 5 i-/ l ^ +5i 

with Aj = - p for i = 0, . . . , k — 1. 



In case of (3 > 7 and /3 Lp > 7 Lp , we prove that there exists some e > such 

that for all e G]0,e ] follows C(y Mi ( QjQ ,))(e) > C{Y ML ( a , p) )(e). From the lower 

bound (25) for C{Y M L[a,a)){ £ ) an d the upper bound (26) for C(Y M L(a,p))( £ ) 
we get the estimate 



C(^AfL(a,Q))(^) - C(XML(a,p))(£) 

> r fe r ^,. M /c 2 , c 2 c 3 y /2 /»r - 

-t^ e ^ T M^rJ 1 - 

| / c 2 c 2 , C3, y/ 2 T ^ h L 2 p - M , -/,,- 

i=0 ^ 

fc-1 

+ E 



c 3 / 1 - 

( (M~ l T) (fc^ - M^/iff 



(l-M^-^Xl-M 2 ^) 
(M^T)^ (/iff +5i - M^liJ - ^; 7+5i + M^- S 'ht a ^ 



V 



+ 



(i-mV-^i-m 1 /) 

-g4^^^) 1/2 M i £ f ? " 

i=0 \ C3 ' L p / 
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£41 ( (^)' /2 A„^« + ^- + * 

i=0 \ \ C 3,^p / 



i=0 



In the following, we make us of M 1 c lt £q**e^ < h La < c^q^e* and 
£ ->■ 0." 



— _L i — - _L I 

c i,pQ 2p£p < < c i,pQ 2p£p , i- e - we have — >■ and ft,£ a — > as 



min{/3-7,0 Lp -7 L?) } 



Multiplying both sides of (27) with K^- e 2 h Lp : and taking into 

account the assumptions Ap > (5 + 7 and 4p > /3 — 7 + 27 ip results in 



1-q 



min{/3-7,0 if> -7 Lf) } 



^ 2 ^ p 2 (C(Y MHa , a) )(e)-C(Y MHa>p) )(e)) 



> 



T 

— 1/2 



1-MT 



/ ,_. /ir^ - - \ 1/2 ^"Ti. 



V 

/ min{/3 ip - 7ij) ,/3-7> N 



40) 1/2 40) f C 3C2,L P . , j 

min{^-7,^ i?) -7 Lp } 

^ p 5 • (28) 



As a result of (28) follows that in the case of j3 — 7 < (3l p — 7£ p there exists 
some e > such that 

C(yML(a,a))(g) > j (2Q) 
C(YML(a,p))( £ ) 

for all £ e]0, £()]■ 111 the case of (3 — 7 = /3l p — 7l p there exists some £0 > 
such that (29) holds for all e e]0,£ ] if c 2 c 3 > (1 — M 2 ^ S ') 2 C2,l p C3,l p and 

(cf) 2 f > (i-m^) 2 (4X) 2 S^- Fina11 ^ = o(£- 2 ) foi- 

lows from (26). 

In case of (3 < 7 and /3 < 2p, we have to compare the dominating terms 
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as e — > 0. Therefore, we get from the lower bound that 



£-7 

C{Y MLM ){e) > ^- e-i-^Tcfl c 2 , ip cj (M V - l)' 
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+ o(e _:4 -") (30) 



and from the upper bound 



-(0) / \ V2 

4,| p (c 2 c 3 c 2iLp J 



- / *(o) *(<>) („_„_ . „_ . V /2 



C 3 C 2 C 3 ^ c 2C2,L p C 3i i p 



1-M^) 2 cf(l-M^) 



+ 4°U,lJ+o(^^) (31) 



c 3 



«, (1 - M*) 



Making use of these two estimates (30) and (31), this results in the estimate 
(20) where (3l p < 7l p because we require that (3l p — = /? — 7 < 0. 

In general, it follows that C(Y ML ^ a ^)(e) = 0[e p J from the 

upper bound (26) for (5 < 7 and any /3 Lp > 0, 7 ip > I. Further, there is 
an asymptotically optimal choice for the parameter q e]0, 1[ such that the 
computational costs are asymptotically minimal. Calculating a lower bound 
for C{YML[a,p)){ e ) an d taking into account the upper bound (31), we get 

C(Y ML{a , p) )(e) = -^—e- 2 - 1 ^ q^C + o^ 2 ^) (32) 

with some constant C > independent of q and e. Now, we have to find some 
q G ]0, 1[ such that 

_ 2 _2^£ g 2 f . -2-2=1. q 2 p , . 

Ce p = mm Ce p 33 

l-q «e]o,i[ 1-g 

for all < e < 1. Solving this minimization problem leads to 

= 7 -/3 
q 7 -/3 + 2p' 

which is asymptotically the optimal choice for q e ]0, 1[ in case of /3 < 7. 
In case of /3 = 7, we get the following lower bound 

c(fM iM )( £ ) > ~t~~ ^ 2 f E $U,o/# + E 4 f V2 L a /* 

1 ~ ^ \i=0 i=0 V C 3:0 / 

+ g (0) / C 2 C 2 , oC 3,o \ V2 ^ + * / C 2 C 2 , C 3 , \ V2 - fr* 

3 V c 3 / a ~t 3 Vc 3 / M Si - 1 



(34) 
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+ c ( 3 'c 2 L a + ^c ( 3 'c 2 L a MSi _° 



1=1 



? v=o V i=o ^ C3,0 ' 

(0) / C2C 2 ,oC 3 ,o \ * (i) T*-fe* ' 



log^- 1 ) log(g^ Cl , a T") \ (o) f log( £ - 1 ) 



alog(M) alog(M) J 3 z \a\og(M) 
,(o) log(£ _1 ) log(g-5c ljQ T Q ) . (0) / log(g-5 Clja T a ^ 



2 



+ ° 3 ° 2 a 2 (log(M))2 + c 3 c 2 ^ a i g(M) 
A / c 2 c 2 c 3 o \ 1/2 T 5i - \ 



where cj^ = c { 3 \ c 2 , La = c 2 , c 3 , La = c 3 , (3 La = P and 7La = 7 for f M L( a ,a). 
Next, we calculate for (5 = 7 the upper bound 



r -2 



C(Y ML{a:P) )(e) <Y^ S 



i=0 i=i V c 3 / 

fc /_ _ _ \ V 2 PL„—tLr, 



i=0 



C3,0 



^ C 3 / i= o ' V C 3,0 / 

+ 6 (o) /c^c^y/^^ 



c 3 



+ E 41 ( ) + E 4°c 2 a, 

i=0 V C 3,^p / i=l 

logte' 1 ) log(g - 5 Ci iP Tp) 



plog(M) plog(M) 



i(0) 



log^- 1 ) , \og(q-^ Cl , p TP) 



2 



^ 3 ° 2 1 plog(M) + plog(M) 

* /~\ /CoCn t Ci r \ X /2 PLy-lLp 



i=0 V \ C 3,L P / 
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:m- l t) s *-i - tfr 1 



i=0 



(36) 



where we applied the relation (24). 



Suppose that (3l p > 7l p and 7, 7l p < 2p. Then, we get from the upper bound 
(36) that C{Yml{<x,p)){^) = 0(e~ 2 (\og(e)) 2 ). Further, comparing the lower and 
the upper bounds (35) and (36), we asymptotically obtain that 



Um C&M)fc) > lim ^- 2 4 0) c 2 (^) Vo(^(log(e)) 2 ) £ 



~° C(Y ML{a , p) )(e) -0 _L_ £ -2 d f) C2 ffl^y + 0(e -2( tog ( e))2) <* 

(37) 

which proves statement (19). □ 

Remark 3.2 Especially, if C3 = C3 and Cz,l v = c^l p , then it follows in case 
of (3 < 7 and (3 < 2p that 

lim C < WrtX £ ) > M 7-/3 f 1 - M V f 1 - ( C2C3 "I V ^ ) 2 . (38) 

TTras, i/ C2C3 < C2,L p C3 i L p ^ follows directly that 

lim C fo L(p * ))(e) > 1 . (39) 
~° C(Y ML{a , p) )(e) 



4 Numerical examples in case of SDEs 



For illustration of the improvement that can be realized with the proposed 
modified multi-level Monte Carlo estimator, we consider the problem of weak 
approximation for stochastic differential equations (SDEs) 

Til 

dX t = a(X t ) dt + J2 V(Xt) dB{ (40) 
3=1 

with initial value X to = x G ~R d driven by m-dimensional Brownian motion. 

In the following, we compare for several numerical examples the root mean- 
square errors (5) versus the computational costs for the multi-level Monte 
Carlo estimator Yml proposed in [5,6,7] and described in Section 2 with the 
proposed modified multi-level Monte Carlo estimator Ym L(a,p) described in 
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Section 3. As a measure for the computational costs, we count the number 
of evaluations of the drift and diffusion functions taking into account the di- 
mension d of the solution process as well as the dimension m of the driving 
Brownian motion. 

In the following, we consider on each level I = 0,1, ... ,L an equidistant dis- 
cretization I ni = {to, . . . ,tr} of [to,^ 1 ] with step size hi — Further, we 
denote by Y n = Y tn the approximation at time t n . In case of the multi-level 
Monte Carlo estimator Yml we apply on each level / = 0,1, ... ,L the Euler- 
Maruyama scheme on the grid 7^ given by Y = x and 

m 

Y n+1 = Y n + a(Y n ) h n + J2 ViXn) I {j) ,n (41) 

where h n = hi and /(,-),„ = B J tn+1 - B\ n for n = 0, 1, 1. The Euler- 

Maruyama scheme converges with order 1/2 in the mean-square sense and 
with order a — 1 in the weak sense to the solution of the considered SDE (40) 
at time T [10]. 

On the other hand, for the modified multi- level Monte Carlo estimator Y ML ( a ^ 
the Euler-Maruyama scheme is applied on levels 0, 1, . . . , L p — 1 whereas on 
level L p a second order weak stochastic Runge-Kutta (SRK) scheme RI6 pro- 
posed in [12] is applied. The SRK scheme RI6 on level L p is defined on the 
grid I hLp by % = x , 

m 

Y n+1 = Y n + \ (a(Y n ) + a(T)) h n + \ ]T (b k (T^) 

k=l 

m 

+ E {¥\Yn) + inr?) + \b k (r^)) / W>B 
k=i 

m 

+ ij2(b k (t^)-b k (t^)) 4h n 

k=l 

m 
k=l 

where h n = h Lp and I (k) ^ n = B k n+1 - B k n for n = 0, 1, . . . , £ - 1 with stages 

m 

T = Y n + a(Y n ) h n + J2 V(Y n ) I {j) , n , 

3=1 

Ti fe) =Y n + a(Y n ) h n ± b k (Y n ) yt, t 



(42) 



Y n ±'£b i (Y n 



(43) 



j'=i 
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Fig. 1. Error vs. computational effort for SDE (45) using f(x) 
f(x) = x 2 (right). 

where I( k ,k),n = l{lf k ),n ~ h n) and 

f I 2^(k),n^(j),n ~ VKJ(k),n) if k < j 

\ \{I{k),nI(j),n + VKiI{j),n) if 3 < k 

based on independent random variables I(k), n with P(I^ >n = -^y/h^) = |. 
Thus, we have a = 1 and p = 2 for the modified multi-level Monte Carlo 
estimator Yml(o,p) in the following. Further, for both schemes the variance 
decays with the same order as the computational costs increase, i. e. (3 = 
Pl p — 7 — 1l p — 1- Then, the optimal order of convergence attained by the 
multi-level Monte Carlo method is 0(s~ 2 (log(e)) 2 ) due to Theorem 2.1. For 
the presented simulations, we denote by MLMC EM the numerical results for 
Yml based on the Euler-Maruyama scheme only and by MLMC SRK the re- 
sults for YML(a,p) based on the combination of the Euler-Maruyama scheme 
and the SRK scheme RI6. 



= x (left) and 



(44) 



As a first example, we consider the scalar linear SDE with d — m — 1 given 
by 

dX t = rX t dt + aX t dB t , X = 0.1 , (45) 

using the parameters r = 1.5 and o = 0.1. We choose T = 1 and apply the 
functionals f(x) = x and fix) = x 2 , see Figure 1. The presented simulations 
are calculated using the prescribed error bounds e = 4 _J for j = 0, 1, . . . , 5. 
In Figure 1 we can see the significantly reduced computational effort for the 
estimator Yml(i,2) (MLMC SRK) compared to the estimator Yml (MLMC 
EM) in case of a linear and a nonlinear functional. 

The second example is a nonlinear scalar SDE with d = m = 1 given by 

dX t = \X t + ^X 2 + 1 dt + ^X 2 + 1 dB t , X = 0. (46) 

We apply the functional f(x) = (log(x + \/x 2 + l)) 3 — 6(log(x + \/x 2 + l)) 2 + 
8 \og(x+\/x 2 + 1). Then, the approximated expectation is given by E(f(X t )) = 
t 3 — 3t 2 + It. Here, the results presented in Figure 2 (left) are calculated for 
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Fig. 2. Error vs. computational effort for the nonlinear SDE (46) (left) and SDE (47) 
(right) with non- commutative noise. 

T = 2 applying the prescribed error bounds e = 4 _ - ? for j — 0, 1, . . . , 6. Here, 
the improved estimator Yml(i,2) performs much better than Yml also for non- 
linear functionals and a nonlinear SDE. 



Finally, we consider a nonlinear multi-dimensional SDE with a d = 4 dimen- 
sional solution process driven by an m = 6 dimensional Brownian motion with 
non-commutative noise: 







/ 243 yl 
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(47) 
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with initial condition Xo = (|, |, 1, |) T . Then, the approximated first mo- 
ment of the solution is given by E{XJ T ) = Xq exp(2T) for i = 1,2,3,4. 
The simulation results calculated at T = 1 for the error bounds e = 4 _J 
for j = 0,1,..., 6 are presented in Figure 2 (right). Again, in the multi- 
dimensional non-commutative noise case the proposed estimator Yml(i,2) needs 
significantly less computational effort compared to the estimator Y ML which 
reveals the theoretical results (19) in Proposition 3.1. 



5 Conclusions 

In this paper we proposed a modification of the multi-level Monte Carlo 
method introduced by M. Giles which combines approximation methods of 
different orders of weak convergence. This modified multi-level Monte Carlo 
method attains the same mean square order of convergence like the originally 
proposed method that is in some sense optimal. However, the newly proposed 
multi-level Monte Carlo estimator can attain significantly reduced computa- 
tional costs. As an example, there is a reduction of costs by a factor (p/a) 2 
for the problem of weak approximation for SDEs driven by Brownian motion 
in case of /3 — 7. This has been approved by some numerical examples for 
the case of p = 2 and a = 1 where four times less calculations are needed 
compared to the standard multi-level Monte Carlo estimator. Here, we want 
to point out that there also exist higher order weak approximation schemes, 
e. g. p = 3 in case of SDEs with additive noise [2], that may further improve 
the benefit of the modified multi-level Monte Carlo estimator. Future research 
will consider the application of this approach to, e.g., more general SDEs like 
SDEs driven by Levy processes [3] or fractional Brownian motion [11] and to 
the numerical solution of SPDEs [13]. Further, the focus will be on numerical 
schemes that feature not only high orders of convergence by also minimized 
constants for the variance estimates. 
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